C.................... CODDN.FOR;1 .......... 16-SEP-1994 14:57:22.60 
      SUBROUTINE PRLN(ID,I,J,JP,IPR,PROD,LOSS,VCON,IVLPS)
C...... this program was written to calculate the production and
C...... losses for the species n(2d),o(1d),no,n(4s). it also calc-
C...... ulates the densities of many minor ions and neutral species
C...... from the assumption of chemical equilibrium. the prod @
C...... loss rates, and concentrations of these species can be
C...... printed in the middle section of the program. in the final
C...... section of the program a calculation of the heating efficiency
C...... of the neutral gas is done.
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL ENTION,PNTXIT,PNTION,ANTION,ANTEX,TOTION,TOTEX
      DIMENSION PROD(6),LOSS(6),RTS(199),R42(VD),VCON(VD),LR(22)
     >  ,FEX(2),VN2O(6),COOL(22)
      COMMON/MION/NOP(VD),N2PLS(VD),O2P(VD)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/ENLOSS/ENL(20,VD)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,IION
C....... for the production rates consult subroutine prdint
      COMMON/PRTER/PUVN2P(VD),PUVO2D(VD),PEO1D(VD),PUVO2P(VD)
     >  ,PUOP2P(VD),UVDISN(VD),PSRNO(VD),PO1DSR(VD),PO1SEL(VD)
     >  ,PN2A3S(VD),PDISOP(VD),PUVIOP(VD),PRN2XS(VD),PRN2AP(VD)
     >  ,PRN2BS(VD),DISNP(VD),DISN4S(VD),DISN2D(VD),DISN2P(VD),PRHEP(VD)
     >  ,HEPLUS(VD),NPLUS(VD),PLYNOP(VD),PHOTN(VD),HE10830(VD)
C----- three commons to hold the EUV, AURoral, and TOTal production rates
C-----  G. GERMANY  11/4/91
      COMMON/EPRDNT/ENTION(3,12,VD),PNTXIT(3,12,VD),PNTION(3,12,VD)
      COMMON/APRDNT/ANTION(3,12,VD),ANTEX(3,12,VD)
      COMMON/PTOTV/TOTION(3,12,VD),TOTEX(3,12,VD)
C  for the IBM since EHT is not used
      EHT=0.0

C*****  obtain reaction rates from subr rats to get their densities
       CALL RATS(J,TE(J),TIZ(J),TN(J),RTS)
C++++++++ multiply o+ + n2 rate by factor for vib n2 ++++
      IF(VCON(J).LE.1) VCON(J)=1.0
      IF(VCON(J).GT.50.0) THEN
         WRITE(6,901)
         WRITE(ID,901)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF
 901  FORMAT(/5X,' ** ERROR: N2* RATE FACTOR EXCESSIVE IN PRLN'/)

      RTS(3)=RTS(3)*VCON(J)
      !.. RTS(99) is used in the RVN2PB file for N2+(v) calculation
      RTS(99)=RTS(42)*RTS(10)

C....... he+ calculation oppenheimer et al
C
      ZNE(J)=NR(1,J)+NR(2,J)+NOP(J)+O2P(J)+N2PLS(J)
C............ o+(2p)
      CALL COP2P(0,J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),OP2P
     > ,PUOP2P(J),0.0D0,HEPLUS(J),N4S(J),NNO(J))
C............ o+(2d)
      CALL COP2D(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),OP2D
     > ,ENTION(1,8,J),PNTION(1,8,J)+ANTION(1,8,J),OP2P,HEPLUS(J),
     >  N4S(J),NNO(J))
C
      NMIN=OP2D+OP2P+NPLUS(J)+HEPLUS(J)
C
      IF(JTER.GT.1.OR.I.EQ.1.OR.I.EQ.-1.OR.I.GT.2) GO TO 12
C.......... newton solution for electron density (zne)
      JITER=0
      ZNE(J)=NR(1,J)+NR(2,J)+NOP(J)+O2P(J)+N2PLS(J)+NMIN
      IF(ZNE(J).LT.10.) ZNE(J)=10.0
C.......THE NEXT THREE LINES SHOULD BE COMMENTED TO KILL VIBN2P............
C....... call vibn2p to get effective n2+ + o->n2 + o+ rate
C...      IF(CHIV(j).LT.1.57) CALL VIBN2P(J,3,0,N2N(J),ON(J),ZNE(J),OP2P
C...     > ,OP2D,RTS,PRN2XS(J),PRN2AP(J),PRN2BS(J),O2N(J),NTOT,ZR(J))
C     IF(CHIV(J).LT.1.57) CALL VIBN2P(J,3,0,N2N(J),ON(J),ZNE(J),OP2P
C    > ,OP2D,RTS,PRN2XS(J),PRN2AP(J),PRN2BS(J),O2N(J),NTOT,ZR(J))
      !.. RTS(99) is used in the RVN2PB file for N2+(v) calculation
      R42(J)=RTS(99)
C
 9    DO 10 ITS=1,2
      IF(JITER.EQ.0) H=ZNE(J)*0.001D0
      IF(ITS.EQ.2) ZNE(J)=ZNE(J)+H
C........... n2+
      CALL CN2PLS(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N2PLS(J),ENTION(3,1,J),ENTION(3,2,J),ENTION(3,3,J)
     > ,PNTION(3,1,J)+ANTION(3,1,J),PNTION(3,2,J)+ANTION(3,2,J)
     > ,PNTION(3,3,J)+ANTION(3,3,J),OP2D,OP2P,HEPLUS(J)
     > ,NPLUS(J),NNO(J),N4S(J))
C
C........... o2+
      CALL CO2P(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PO2P,
     > O2P(J),ENTION(2,7,J),PNTION(2,7,J),NR(1,J),OP2D,N2PLS(J),
     > NPLUS(J),NV(2,J),NV(3,J),OP2P)
C
C........... no+
      CALL CNOP(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PNOP,NOP(J)
     >  ,NR(1,J),N2PLS(J),O2P(J),NV(2,J),NV(3,J),NPLUS(J),N2P(J)
     >  ,PLYNOP(J),VCON(J),N2D(J),OP2D)
C
C.....calculation of ne from the total ion concentration
      B=NR(1,J)+NR(2,J)+N2PLS(J)+NMIN
      C=PNOP/RTS(5)+PO2P/RTS(6)
      IF(JITER.EQ.0) ZNE(J)=(B+DSQRT(B*B+4.0*C))/2.0
      FEX(ITS)=ZNE(J)-NR(1,J)-NR(2,J)-NOP(J)-O2P(J)-N2PLS(J)-NMIN
 10   CONTINUE
C
      JITER=JITER+1
      IF(JITER.EQ.1) GO TO 9
      DEX=(FEX(2)-FEX(1))/H
      ZNE(J)=ZNE(J)-H-FEX(1)/DEX
      IF(JITER.GT.9) GO TO 14
      IF(ABS(FEX(1)/DEX/ZNE(J)).GT.1.0E-2) GO TO 9
      DX=ABS(FEX(1)/DEX)
      IF(JITER.GT.5) WRITE(6,58) J,JITER,ZNE(J),DX
 58   FORMAT(3I6,1P,22E10.2)
 14   IF(JITER.GT.22) THEN
         WRITE(6,57) J
 57      FORMAT('   NEWTON NON-CONVERGENCE IN SUB PRLN AT PT. J=',I3)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF     
 12   CONTINUE
C....... Store N2+(v) + O reaction rate for N2+ subroutine
      RTS(99)=R42(J)
C
C........n2(a3sigma)
      CALL CN2A(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N2A,TOTEX(3,1,J),TOTEX(3,2,J),TOTEX(3,3,J),TOTEX(3,4,J))
C...... densities for O2 atmospheric bands
        CALL CO2(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,BOD3,O1D(J),O2SS,TN(J),O2B1,O2DEL,O2ADEL,O2ASIG)
C.......O(1S)
      CALL CO1S(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,O1S,O2P(J),N4S(J),PO1SEL(J),NPLUS(J),N2A,PO1DSR(J),O2ADEL,O2DEL)

      !......... production rates for metastables and odd nitrogen ........

       CALL CN2P(J,0,0,ZR(J),RTS,ON(J),O2N(J)
     > ,N2N(J),ZNE(J),PROD(5),LOSS(5),N2P(J),DISN2P(J),UVDISN(J),O2P(J)
     > ,NV(3,J),N2PLS(J))
      IF(IPR.LT.5) N2P(J)=PROD(5)/LOSS(5)
C
      CALL CN2D(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,NOP(J),ZNE(J),PROD(1),LOSS(1),N2PLS(J),DISN2D(J),UVDISN(J)
     > ,NPLUS(J),N2P(J),N2D(J),NR(1,J),NNO(J),N2A)
C
      IF(IVLPS.EQ.-9) N2D(J)=PROD(1)/LOSS(1)
      IF(IVLPS.EQ.-9) NV(1,J)=PROD(1)/LOSS(1)
C
      CALL CN4S(J,0,0,ZR(J),RTS,ON(J)
     > ,O2N(J),N2N(J),ZNE(J),PROD(2),LOSS(2),N4S(J),DISN4S(J),NV(1,J)
     > ,N2P(J),NR(1,J),N2PLS(J),UVDISN(J),NOP(J),NPLUS(J),NV(3,J)
     > ,O2P(J),PSRNO(J),VCON(J))
C
      CALL CNO(J,0,0,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     >   ,ZNE(J),PROD(3),LOSS(3),NV(1,J),NV(2,J),N2P(J)
     > ,NNO(J),O2P(J),NR(1,J),PSRNO(J),PLYNOP(J),N2A,NPLUS(J))
C
      CALL CO1D(0,J,0,0,ZR(J),RTS,ON(J),O2N(J)
     > ,N2N(J),ZNE(J),PROD(4),LOSS(4),O2P(J),NV(1,J),PEO1D(J),PO1DSR(J)
     > ,O1S,NPLUS(J),O1D(J))
C
      IF(JTER.LE.1.AND.IPR.LT.3) NNO(J)=PROD(3)/LOSS(3)
      IF(JTER.LE.1.AND.IPR.LT.4) O1D(J)=PROD(4)/LOSS(4)
      IF(I.LT.2) RETURN
C
C
C...... this next section for printing densities and production rates
C
      IF(IVLPS.EQ.-9) O1D(J)=PROD(4)/LOSS(4)
C
      IF(I.EQ.2) CALL CION(J,ID,JP,ZR(J),NR(1,J),O2P(J),NOP(J),N2PLS(J)
     > ,NPLUS(J),OP2D,OP2P,ZNE(J),TE(J),TIZ(J),VCON(J),HEPLUS(J)
     > ,NR(2,J),NNO(J),TN(J))
C
      IF(I.EQ.3) CALL CNEUT(J,ID,JP,ZR(J),O1S,N2A,N2D(J),O1D(J),NNO(J),
     > N4S(J),N2P(J))
C
C------ Some production rates for printing  converts SP to DP
      P3X11=TOTEX(3,11,J)
      P3X12=TOTEX(3,12,J)
      P3X2=TOTEX(3,2,J)
      P3X3=TOTEX(3,3,J)
      T3I6=TOTION(3,6,J)
C
C----- important emission rates
      IF(I.EQ.4) CALL CEMM(J,ID,JP,ZR(J),RTS,O1D(J),N2D(J),OP2D,OP2P,O1S
     > ,DISNP(J),ZNE(J),ON(J),N2N(J),N4S(J),P3X12,P3X11,N2A,P3X2,P3X3
     > ,T3I6)
C
      IF(I.EQ.5) CALL CO2EM(J,ID,JP,ZR(J),O2B1,O2DEL,O2SS,O2ADEL,O2ASIG)
C
      IF(I.EQ.6) CALL CN2D(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),NOP(J)
     > ,ZNE(J),PP,LL,N2PLS(J),DISN2D(J),UVDISN(J),NPLUS(J),N2P(J)
     > ,N2D(J),NR(1,J),NNO(J),N2A)
C
C ------ O(1D) production rates
      IF(I.EQ.7) CALL CO1D(0,J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),PP,LL,O2P(J),N2D(J),PEO1D(J),PO1DSR(J),O1S,NPLUS(J)
     > ,O1D(J))
C
      !.. 6300 Volume emission rates
      IF(I.EQ.36) CALL CO1D(1,J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),PP,LL,O2P(J),N2D(J),PEO1D(J),PO1DSR(J),O1S,NPLUS(J)
     > ,O1D(J))
C
      IF(I.EQ.8) CALL CNO(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PP,LL,N2D(J),N4S(J),N2P(J),NNO(J),O2P(J),NR(1,J),PSRNO(J)
     > ,PLYNOP(J),N2A,NPLUS(J))
C
      IF(I.EQ.9) CALL CN4S(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PP,LL,N4S(J),DISN4S(J),N2D(J),N2P(J),NR(1,J),N2PLS(J)
     > ,UVDISN(J),NOP(J),NPLUS(J),NNO(J),O2P(J),PSRNO(J),VCON(J))
C
      IF(I.EQ.10) CALL CO1S(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,O1S,O2P(J),N4S(J),PO1SEL(J),NPLUS(J),N2A,PO1DSR(J),O2ADEL,O2DEL)
C
      IF(I.EQ.11) CALL CN2PLS(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),N2PLS(J),ENTION(3,1,J),ENTION(3,2,J),ENTION(3,3,J)
     > ,PNTION(3,1,J)+ANTION(3,1,J),PNTION(3,2,J)+ANTION(3,2,J)
     > ,PNTION(3,3,J)+ANTION(3,3,J),OP2D,OP2P,HEPLUS(J)
     > ,NPLUS(J),NNO(J),N4S(J))
C
      IF(I.EQ.12) CALL CNOP(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PNOP,NOP(J),NR(1,J),N2PLS(J),O2P(J),N4S(J),NNO(J),NPLUS(J)
     > ,N2P(J),PLYNOP(J),VCON(J),N2D(J),OP2D)
C
      IF(I.EQ.13) CALL CO2P(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PO2P,O2P(J),ENTION(2,7,J),PNTION(2,7,J),NR(1,J),OP2D,N2PLS(J),
     > NPLUS(J),N4S(J),NNO(J),OP2P)
C
C---- Calculate O+(4S) production and loss rates
      IF(I.EQ.14) THEN
         EUVP4S=ENTION(1,7,J)
         PEOP4S=PNTION(1,7,J)+ANTION(1,7,J)
         CALL COP4S(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),NR(1,J)
     >    ,EUVP4S,OP2D,OP2P,PEOP4S,PDISOP(J),N2PLS(J),N2D(J),NNO(J)
     >    ,VCON(J),HEPLUS(J),NPLUS(J))
      ENDIF
C
      IF(I.EQ.15) CALL COP2D(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),OP2D,ENTION(1,8,J),PNTION(1,8,J)+ANTION(1,8,J),OP2P,
     >  HEPLUS(J),N4S(J),NNO(J))
C
      IF(I.EQ.16) THEN
         PEOP2P = PNTION(1,9,J)+ANTION(1,9,J)
         CALL COP2P(1,J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),OP2P
     >    ,PUOP2P(J),PEOP2P,HEPLUS(J),N4S(J),NNO(J))
      ENDIF
      !.. 7320 Volume emission rates
      IF(I.EQ.22) THEN
         PEOP2P = PNTION(1,9,J)+ANTION(1,9,J)
         CALL COP2P(1,J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),OP2P
     >    ,PUOP2P(J),PEOP2P,HEPLUS(J),N4S(J),NNO(J))
      ENDIF

      PEDISNP= PNTION(3,4,J)+PNTION(3,5,J)+PNTION(3,6,J)
      IF(I.EQ.17) CALL CNPLS(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),DISNP(J),PEDISNP,NPLUS(J),NR(1,J),N2D(J),OP2P,HEPLUS(J)
     > ,PHOTN(J),O2P(J),N4S(J),OP2D,N2PLS(J),NNO(J))
C
      IF(I.EQ.18) CALL CN2A(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N2A,TOTEX(3,1,J),TOTEX(3,2,J),TOTEX(3,3,J),TOTEX(3,4,J))
C
      IF(I.EQ.19) CALL CN2P(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PP,LL,N2P(J),DISN2P(J),UVDISN(J),O2P(J),NNO(J),N2PLS(J))
C
      IF(I.EQ.20) CALL CHEP(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,HEPLUS(J),PRHEP(J),HE(J),NNO(J),NR(1,J),NR(2,J),HN(J))
C
      IF(I.EQ.21) CALL GRAV_WAV(J,ID,JP,ZR(J),ON(J),O2N(J),N2N(J),HE(J)
     > ,HN(J),TN(J))
C
      !.. IF(I.EQ.22) CALL CRT2(J,ID,JP,ZR(J),RTS) !.. reaction rates not used
C
      IF(I.EQ.23) CALL CMSIS(J,ID,JP,ZR(J),ON(J),O2N(J),N2N(J),HE(J)
     > ,HN(J),TN(J),O1S,N2A,N2D(J),O1D(J),NNO(J),N4S(J),N2P(J))
C
      IF(I.EQ.24) CALL C1200(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),P3X12,P3X11,N4S(J),DISNP(J))
C
      IF(I.EQ.25) CALL CO2(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,BOD3,O1D(J),O2SS,TN(J),O2B1,O2DEL,O2ADEL,O2ASIG)
C
      IF(I.EQ.26) CALL CN2PV(J,ID,JP,ZR(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),OP2P,OP2D,PRN2XS(J),PRN2AP(J),PRN2BS(J))
C
      IF(I.EQ.27) CALL PRDPRA(J,ID,JP,ZR(J))
C
      !--------------------CGAG  G. GERMANY  11/22/91---------------------|
      !.CGAC.ADDED CALL TO CEMOX USING INTERPOLATED VALUES

      IF(I.EQ.28) THEN

      !.CGAC.FIRST CALCULATE PASSED VALUES: SUM OF AURORAL AND PE SOURCES
       P1=TOTEX(1,1,J)
       P2=TOTEX(1,2,J)
       P3=TOTEX(1,3,J)
       P4=TOTEX(1,4,J)
       P5=TOTEX(1,5,J)
       P6=TOTEX(1,6,J)
       P7=TOTEX(1,7,J)
       P8=TOTEX(1,8,J)
       P9=TOTEX(1,9,J)
       P10=HE10830(J)
      !.. Electron impact excitation rates for atomic oxygen.Modified 2020-08-02 
      !.. to include O1D to calculate 6300 emission rates
       CALL CEMOX(J,ID,JP,ZR(J),RTS,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,
     >   ZNE(J),NR(1,J),O1D(J),O1S)

      ENDIF

      !****************************************************
      !.. the following section calculates the various components
      !.. of the euv heating efficiency of the neutral atmosphere
      !*******************************************************

      !... Electron cooling
 45   IF(I.NE.29) GO TO 46
      !.. NMINOR=NOP(J)+O2P(J)+N2PLS(J)
      !..CALL ECRATS(ID,J,ZR(J),TE(J),TIZ(J),TN(J),NR(1,J),NR(2,J)
      !..> ,NMINOR,ON(J),N2N(J),O2N(J),HN(J),HE(J),EHT,TLSS,FOL,COOL)
      !..ENL(8,J)=COOL(5)
      !..ENL(20,J)=FOL
      IF(JP.EQ.1) WRITE(ID,191) 
 191  FORMAT(/'  ** Electron heating and cooling not available on'
     >  ' this file **'/)
 46   CONTINUE

      IF(I.NE.30) GO TO 51
      IF(JP.EQ.1) WRITE(ID,150)
 150  FORMAT(/3X,'FOR HEATING EFFICIENCY: O2 dissociative energy'
     > ,1X,'losses '
     > /,8X,'---------- Non 3-body dissociation ---------- '
     > '  --- 3 body not in heating effic ----'
     > ,/3X,'ALT   S-R    N2D+O2  N4S+O2  N++O2   E+O2+'
     > ,3X,'N+O2+    O2**   O2del    O2b1   O2ADEL  O2ASIG  total'
     >  ' total3body')
      LR(1)=5.08*PO1DSR(J)
      LR(2)=5.08*N2D(J)*O2N(J)*RTS(16)
      LR(3)=5.08*RTS(7)*O2N(J)*N4S(J)
      LR(4)=5.08*RTS(25)*O2N(J)*NPLUS(J)
      LR(5)=5.08*RTS(6)*ZNE(J)*O2P(J)
      LR(6)=5.08*RTS(21)*O2P(J)*N4S(J)
      LR(7)=2.3*O2SS*(5.0E-13*O2N(J)+3.0E-11*ON(J))
      LR(8)=O2DEL*(1.3E-16*ON(J)+2.6E-20*TN(J)**.78*O2N(J)+1.0E-20*
     > N2N(J))
      LR(9)=.7*O2B1*2.2E-15*N2N(J)
      LR(10)=O2ADEL*3.0*2.0E-13*N2N(J)
      LR(11)=3.0*O2ASIG*(3.0E-13*N2N(J)+3.0E-13*ON(J))
      ENL(1,J)=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)
      ENL(4,J)=LR(7)+LR(8)+LR(9)+LR(10)+LR(11)
      WRITE(ID,88) ZR(J),(LR(K),K=1,11),ENL(1,J),ENL(4,J)
      RETURN
C
 51   IF(I.NE.31) GO TO 52
      IF(JP.EQ.1) WRITE(ID,151)
 151  FORMAT(/2X,'FOR HEATING EFFICIENCY: Forbidden radiation losses'
     > /1X,'ALT  O(1D)   O(1S)  O+(2P)   N(2A)   N(2P)   N(2D)   O2atm'
     > ,3X,'N1200   total')
      LR(1)=1.96*O1D(J)*RTS(54)
      LR(2)=2.05*(RTS(55)+RTS(56))*O1S
      LR(3)=1.70*0.171*OP2P
C...... N(2A) emissions result in 4 not 6eV due to vibrational excitation
      LR(4)=4.0*0.57*N2A
      LR(5)=(3.58*RTS(58)+1.119*RTS(57))*N2P(J)
      LR(6)=2.4*N2D(J)*RTS(61)
C...... O2 atmospheric O+O+M
      LR(7)=1.7*0.083*O2B1
      LR(8)=10.33*P3X11*N4S(J)
      ENL(2,J)=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)+LR(8)
      WRITE(ID,88) ZR(J),(LR(K),K=1,8),ENL(2,J)
      RETURN
C
 52   IF(I.NE.32) GO TO 53
      IF(JP.EQ.1) WRITE(ID,152)
 152  FORMAT(/2X,'FOR HEATING EFFICIENCY: Permitted radiation losses. '
     > ,2X,'63u(e*) 63u(O+O) 5.3u(NO) are not included in heating'
     > ,1X,'efficiency calculation'
     > ,/1X,'ALT  O+(4P) O+(2P*)   1pos    2pos    LBH     O2+     N2+'
     > ,1X,'     O2     O    63u(e*) 63u(O+O) 5.3u(NO)   total')
C,,,,,,,, for atomic oxygen ions >o+(4p)*
       LR(1)=11.71*(TOTION(1,4,J))
       LR(2)=19.20*(TOTION(1,5,J))
C.......... heat radiated by n2(b->a, 0.5c->b->a, a->x, others)
       LR(3)=1.1*TOTEX(3,2,J)
       LR(4)=3.9*TOTEX(3,3,J)
       LR(5)=2.3*TOTEX(3,4,J)
C........... o2+ (1st 4 states)
      LR(6)=4.0*TOTION(2,2,J)+ 4.8*TOTION(2,3,J)+ 6.1*TOTION(2,4,J)
C........ radiation from n2+ ......
      LR(7)=1.15*TOTION(3,2,J)+ 3.2*TOTION(3,3,J)
C,,,,,,,,,,, o2(1st four levels)
      LR(8)=5.0*TOTEX(2,12,J)
C...... radiations from O minus O(1D) and O(1S)
      LR(9)=12.0*TOTEX(1,12,J)
C........ e* + O -> 63u
      LR(10)=ENL(20,J)
      ENL(3,J)=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)+LR(8)+LR(9)
     > +LR(10)
C
C,,,,,,,,,,,, 63 micron radiation ---- o-o collisions
      ENL(6,J)=ON(J)*1.044E-6*EXP(-228.0/TN(J))/
     >       (1.0+0.6*EXP(-228.0/TN(J))+0.2*EXP(-325.0/TN(J)))
C
C........ cooling by 5.3 micron (no*) kockarts grl february 1980 p137
      W=6.5E-11*ON(J)/(6.5E-11*ON(J)+13.3)
      ENL(13,J)=0.232875*NNO(J)*W*13.3*EXP(-2700.0/TN(J))
      LR(11)=ENL(6,J)
      LR(12)=ENL(13,J)
      WRITE(ID,88) ZR(J),(LR(K),K=1,12),ENL(3,J)
      RETURN
 53   IF(I.NE.33) GO TO 55
      IF(JP.EQ.1) WRITE(ID,154)
 154  FORMAT(/2X,'FOR HEATING EFFICIENCY: Metastable kinetic heating'
     > ,/1X,'ALT  N2D+O  N2D+O2   N2D+e  O1D+N2  O1D+O2  O+2D+N2  O+N2+'
     > ,1X,' e+O+2D  O+O+2P'
     >,'  N2+O+2P e+O+2P   O+N2A   O+N2P   O+N2v   O2 qch   total')
      LR(1)=2.38*N2D(J)*ON(J)*RTS(15)
      LR(2)=1.84*N2D(J)*O2N(J)*RTS(16)
      LR(3)=2.38*N2D(J)*ZNE(J)*RTS(8)
      LR(4)=1.96*O1D(J)*N2N(J)*RTS(33)
      LR(5)=1.96*O1D(J)*O2N(J)*RTS(34)
      LR(6)=1.33*OP2D*RTS(19)*N2N(J)
      LR(7)=2.3*N2PLS(J)*ON(J)*R42(J)
      LR(8)=3.31*OP2D*ZNE(J)*RTS(12)
      LR(9)=5.00*RTS(26)*ON(J)*OP2P
      LR(10)=3.02*RTS(20)*N2N(J)*OP2P
      LR(11)=1.69*RTS(13)*ZNE(J)*OP2P
      LR(12)=6.14*RTS(36)*ON(J)*N2A
      LR(13)=3.0*RTS(37)*N2P(J)*ON(J)
C....... n2 vibrational heating via o and o2 quenching
      LR(14)=0.0D0
      EXPTN=EXP(3380.0/TN(J))
      C01=TN(J)**2.87*2.4685E-22/EXPTN
      C01ON=C01*ON(J)
      DO 379 IS=2,6
C******* B7=prod due to t-v between o and n2
      IF(IS.NE.1) B7=C01ON*((IS-1)*STR(IS-1,J)+IS*EXPTN*STR(IS+1,J))
C******* B8=loss due to quenching of n2 by o
      IF(IS.LE.9) B8=C01ON*((IS-1)*EXPTN+IS)*STR(IS,J)
      LR(14)=LR(14)+0.28*(B8-B7)
 379  CONTINUE
       IF(LR(14).LT.0) LR(14)=0.0D0
      LR(15)=ENL(4,J)
      ENL(4,J)=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)+LR(8)
     > +LR(9)+LR(10)+LR(11)+LR(12)+LR(13)+LR(14)
      WRITE(ID,88) ZR(J),(LR(K),K=1,15),ENL(4,J)
      RETURN
C
 55   IF(I.NE.34) GO TO 56
      IF(JP.EQ.1) WRITE(ID,155)
 155   FORMAT(/2X,'FOR HEATING EFFICIENCY: Ground state kinetic heating'
     > ,/1X,'ALT  O+N2+   E+N2+   O2+N2+  E+NO+   E+O2+   N+O2+   N+O2+'
     > ,1X,'  NO+O2+  N2+O+   O2+O+   N+NO    N+O2    O2+N+   O2+N+'
     > ,1X,'   O+N+    total')
      LR(1)=0.70*RTS(10)*ON(J)*N2PLS(J)
      LR(2)=3.44*RTS(11)*ZNE(J)*N2PLS(J)
      LR(3)=3.53*RTS(17)*O2N(J)*N2PLS(J)
      LR(4)=0.90*NOP(J)*ZNE(J)*RTS(5)
      LR(5)=4.42*RTS(6)*ZNE(J)*O2P(J)
      LR(6)=4.19*RTS(21)*O2P(J)*N4S(J)
      LR(7)=0.05*RTS(21)*O2P(J)*N4S(J)
      LR(8)=2.81*RTS(23)*O2P(J)*NNO(J)
      LR(9)=1.10*NR(1,J)*N2N(J)*VCON(J)*RTS(3)
      LR(10)=1.55*NR(1,J)*O2N(J)*RTS(4)
      LR(11)=1.85*RTS(9)*N4S(J)*NNO(J)
      LR(12)=1.42*RTS(7)*O2N(J)*N4S(J)
      LR(13)=6.67*RTS(30)*O2N(J)*NPLUS(J)
      LR(14)=0.11*RTS(25)*O2N(J)*NPLUS(J)
      LR(15)=0.93*RTS(31)*ON(J)*NPLUS(J)
      ENL(5,J)=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)+LR(8)
     > +LR(9)+LR(10)+LR(11)+LR(12)+LR(13)+LR(14)+LR(15)
      WRITE(ID,88) ZR(J),(LR(K),K=1,15),ENL(5,J)
 56    CONTINUE
C........ Heating from 3 body collisions O+O+M
      ENL(7,J)=BOD3*5.08
C
C........ ENL(8,J) = thermal e cooling. See sub ECRATS on READMA.FOR

C........o2 dissoc heating. Energy left after producing O(1D)
      ENL(9,J)=ENL(11,J)-7.04*PO1DSR(J)
      IF(ENL(9,J).LT.0.0) ENL(9,J)=0.0
C
C........ initial energy partition (euv deposition rate) - first into ions
       ENL(10,J)=ENTION(1,10,J)
C........ energy to photoelectrons (and auroral electrons)
       ENL(15,J)=ENTION(1,11,J)+ANTION(1,11,J)
C
C........ ENL(11,J)= S-R deposition rate. See sub SCHUMN on RSPRIM.FOR
C
C.............. kinetic heating from n2 euv dissociation.........
       ENL(12,J)=1.0*UVDISN(J)+2.0*DISNP(J)

C....If I = 40, 41, 42 or 43 calculate hot O population.
        IF (I.GE.40.AND.I.LE.44) THEN

C------ Hot oxygen production rates from N2*
          EXPTN=EXP(3380/TN(J))
          C01=TN(J)**2.87*2.4685E-22/EXPTN
          C01ON=C01*ON(J)
          DO 656 IL=1,6
            VN2O(IL)=C01ON*((IL-1)*EXPTN+IL)*NVS(IL,J)
 656      CONTINUE

          CALL HOTOXY(ID,J,I,JP,ZR(J),RTS,NOP,O2P,ON,N4S,N2N,O2N,ZNE
     >    ,N2P,N2D,O1D,NNO,NPLUS,OP2D,OP2P,NR,VN2O,HN,TIZ,HEPLUS(J))
        END IF
C
 88   FORMAT(F6.1,1P,22E8.1)
 89   FORMAT(I4,1P,8E8.1,0PF6.0,F4.1)

      RETURN
       END
